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Abstract 

We study a reaction model that presents stochastic resonance purely due to internal noise. This 
means that the only source of fluctuations comes from the discrete character of the reactants, and 
no more noises enter into the system. Our analysis reveals that the phenomenon is highly complex, 
and that is generated by the interplay of different stochasticity at the three fixed points of a bistable 
system. 
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Stochastic resonance is a phenomenon that has been extensively studied in the literature, 
both for its inherent interest, as well as for its broad range of applications [1]. Since the 
seminal works on this subject appeared Q], much theoretical work has been done in order to 
extend the variety of situations in which similar phenomena arise. Very interestingexam- 
ples, as stochastic coherence 0], the effect of coloured noise multiplicative noise p, and 
noise-mediated localization |6[ have already been found and analyzed. A different and very 



relevant question is how to get some type of stochastic resonance purelv due to internal fluc- 
tuations. One answer was given in the form of system size resonance 7(, obtained when the 
optimal output of a system is achieved for a certain finite number of constituent subsystems. 
Stochastic optimization also appears by tuning a continuous parameter: numerical studies 
have shown that chemical reactions can undergo stochastic coherence due to internal noise, 

n q 

both in autonomous [8] and non-autonomous situations [9|. But however it is remarkable the 
seemingly little attention that has been paid on internal noise amplification of an external 
signal in chemical kinetics mimicking the classical situation in stochastic resonance, apart 
from the seminal work by Dykman et al. . The main goal of the present work is to study 
this phenomenon by means of analytical tools, and without relying on the assumption of 
Gaussian noise, what allows us to study the problem even if one state is almost empty. 
Reaction kinetics is a prototypical problem where the effect of the fluctuations can be 



quantified Internal fluctuations in these systems appear due to the discrete character 
of the reactants, and a mean-field description always omits some of their features. We 
will show that stochastic resonance is one of them. For our purposes we will consider the 
following series of reactions: — > X, at rate K\, X — > 0, at rate K 2 , 2X —>■ 3X, at rate 
K 3 , and 3X — > 2X, at rate K4. We can interpret this problem as a population dynamics 
model: the first reaction represents emigration into the system, the second death, the third 
reproduction, and the last competition for limited resources of nutrient. A similar problem 
but without the three body reaction (the competition reaction) was studied in Ref. jl^ . 
and we will show here that including this new reaction opens the possibility of stochastic 
resonance. This system can be modelled by means of a combinatorial master equation 

= Kl [ P ( n _ M ) _ p( n> t )] + K 2 [(n + l)P(n + 1, t) - nP(n, t)} + 
-y [{n - l)(n - 2)P(n -l,t)-n(n- l)P(n, £)] + 
-y [(n + l)n(n - l)P(n + 1) - n(n - l)(n - 2)P(n, t)}. (1) 



To study this master equation we will employ the techniques developed by Elgart and 



Kamenev 



121] . Consider the generating function 

oo 

G(p,t) = J2p n P(n,t), (2) 

n=0 

this function obeys the imaginary time Schrodinger equation d t G = —HG, where the Hamil- 
tonian is given by 

H(p q) = K 1 (l-p) + K 2 (p - l)q + ^(1 - p)p 2 q 2 + ^ (p - l)p 2 f. (3) 
The effect of the momentum and coordinate operators is pG(p,t) = pG(p,t) and qG(p,t) = 



—d p G(p,t). The detailed procedure is described in |12j], but it is worth pointing out that 
similar techniques were previously introduced in the literature Q]. We can write now the 
classical equation of motion, that is given by 



1 = ~n 



dq dH 



= *i - K 2 q + ^q 2 - ^g 3 , (4) 
P =i 2 6 



dt dp 

in which the coordinate q plays the role of mean-field density. In order to understand the 
phenomenon of stochastic resonance in a chemical system we need some kind of bistability. 
For this, we make the next assumptions on the rate constants: K\ = eLn 3 , K 2 = n 2 (e + 
L + eL), K3 = 2n(l + e + L), and K4 = 6, where 0<e<l,n>0, and L > 1. This way 
Eq. (JU) has three fixed points: g_ = en, q = n, and q + = Ln, q_ and q + are stable and 
go is unstable. These new parameters e and L are auxiliary mathematical variables that 
are introduced in order to simplify the notation; they have, in contrast, no direct physical 
meaning but the one given through their relation with the rate constants. The variable n 
is a measure of the size of the system, and will be the tuning parameter that will allow 
us to find the resonance, because we can modify its value without changing the relative 
distance between the fixed points. On the other hand, assuming bistability is quite natural 
in population dynamical systems, as for instance in those modelling epidemics Q|. 

The deterministic equation predicts an evolution to one of the stable fixed points, g_ 
for initial conditions g(0) < q , and q + if g(0) > q . In the stochastic case, both states 
become metastable, and the system will jump from one to the other indefinitely. In this 
case the system does not experience the same noise strength in the neighbourhood of the 
fixed points, since the noise has not been imposed externally in an ad hoc manner. To 
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estimate the frequency of jump we will employ the instanton technique for reaction kinetics, 



developed in the paper by Elgart and Kamenev |12| . Let us briefly explain why we do think 
that this method is the most appropriate for the present case. It is known that this type of 
theories, as the one given by Hamiltonian (j3J), can be mapped onto a stochastic differential 
equation jl5|; so one could be suggested to derive the stochastic equation and then apply 
standard techniques [l]. However, the presence of cubic powers of the momentum in the 
Hamiltonian denote the presence of non-Gaussian noise, and its exact correspondence with 
Langevin equations is still unknown 

The action corresponding to Hamiltonian Eq. (JHJ) is 

S[p,q]=Et- [ qpdt-q(0)\p(0)-l], (5) 
Jo 

where we have used the fact that E = H(p, q) is an integral of motion, i. e. E — 0. To 
simplify our calculations we will perform the change of variables q = qn and t = t/n 2 , so 
the action becomes 

S\p, q] = nH(p, q)t — n / qpdt — nq(0)[p(0) — 1], (6) 

Jo 

where H denotes the Hamiltonian with the new rate constants, and the fixed points have 
now the values g_ = e, q + = L, and go = 1. The tildes will be suppressed from now on. 

In the long time limit, t — > oo, the system approaches the trajectories of zero energy, 
H(p,q) = 0. The decay time from one state to the other is estimated as r = aexp(5'o), 
where a is the relaxation time to the arrival state and So is the action along the non-mean- 
field line [3]. Let us briefly comment on this point. The equation H(p,q) = has two 
solutions, p = 1, the mean-field line, and the non-mean-field or instanton line, both are 
depicted in Fig. Qfor e = 1/2, n = 1, and L = 3. One can see that if we only consider the 
mean-field line, then q_ and q + are stable and go is unstable. But if we consider the whole 
phase space, then we realize that the three points are actually saddles. The time it takes 
to go from one point to the other is proportional to the exponential of the action along the 
line, which is identically zero for mean-field lines. Since we are performing a semiclassical 
calculation, it is valid only for calculating transition times along non-mean-field lines if 
Sq > 1. We will estimate the transition time from g_ to g + (and vice versa) as the time it 
takes to move along the non-mean-field line, since this time is much longer than the time it 
takes to perform the mean-field part of the trajectory. 
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FIG. 1: Phase space for e = 1/2, n = 1, and L = 3. The mean-field line is shown dashed and the 
non-mean-field line is represented as a solid line. The fixed points are encircled and the arrows 
show the direction of motion along the different lines. Note that if we only consider the mean-field 
line, then q = 1/2 and q = 3 are stable and q = 1 is unstable. However, taking into account the 
whole phase space the three fixed points become saddles. 

In transitions from g_ to go the action is So = n(g_ — qo) + n pdq, where p can be 
obtained from the relation H(p, q) = 0: 



Although this integral can be computed straightforwardly, the resulting expression is cum- 
bersome, and would not help us to understand the phenomenon. So we will study the simpler 
case e — > 0: 



The reader might wonder how the system can leave the empty state since there is no emi- 
gration. This is true, because in an empty system there are no fluctuations, and once there, 
we will stay in it forever. However, the empty state might be metastable (it is very easy to 
find a set of reactions that rends the empty state unstable), and a perturbation could lead 
us out of it. So in this case e should be considered very small (e << 1), but not identically 




(7) 



limS'o = n(g_ — go) + 2ny/L ■ arccsc ( Vl + L ) = 2n\/~L 





zero. Now it appears clear another advantage of the Elgart-Kamenev approach: it allows us 
to estimate the effect of internal fluctuations in a state with almost no population, contrary 
to perturbative approaches [llj]. A simpler expression can be obtained by taking the limit 
L — > oo: 5*o = n(q_ — go) + 2n = n. Note that in this case the action reduces to the distance 
between the fixed points, suggesting that this distance is the key parameter that rules the 
transitions between the metastable states. We will see, however, that the dynamics is not 
as simple as this. 

We can also compute the action along the non-mean-field line from q + to go 
So = n(q + — qo) + n / pdq — > n(L — 1) + 2n\/L arctan ( yf~L J — arccot ( V~L ) , e — ► 0. 

Jq+ 

(9) 

In the limit L — > oo the second term in the right hand side goes like ~ y/L, so it is irrelevant 
compared to the first term, that is the difference between g and q + . So we have arrived at 
the same conclusion as in the last paragraph. 

Now, if we want to spend the same time in both trips g_ — ► go and q + — ► go we need to 
identify expressions Eq.(JEJ) and Eq.(jSJl, so we are led to solve the transcendental equation 
y/Z/2 + arccot (v^j = arccsc (y/L + l) + arctan (y/L^. This equation can be solved nu- 
merically to get L* m 5.43. At this point one is tempted to use the simplified expressions 
obtained in the limit L — > oo in order to get a closed expression for L*. However, equalling 
the most relevant terms in this limit, one arrives at the contradiction L* — 1. This means 
that we have to take into account that the three points are at a finite distance from each 
other, because the transition times depend on the relative position of the three fixed points 
altogether. Note that assuming the equality of both decay times is not essential for the 
system to undergo resonant behaviour, but it simplifies the forthcoming calculations. 

So we finally have our bistable system (g_,go,g+) = (0, 1,L*). This system is extremely 
asymmetric, not only because the distance between the fixed points is not same, but more 
importantly because the "potential energy" at these points is very different. We have defined 
the potential energy as the integral of the force V(q) — — J F(q)dq, where the force is given 
by F(q) = —d g H(p,q)\ p= i. We have depicted the potential in Fig. El where one can see 
the huge difference of depth between the two wells. From this figure one realizes that the 
fluctuations are much stronger in the case of a higher population, and so, they have to 
be stopped by a higher energy barrier. The fluctuations are so weak in the case of a low 
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FIG. 2: Potential energy in the system with e = 0, n = 1, and L = L* as explained in the text. 
The minima are displayed at q = and q = L*, and the maximum is located at q = 1. One can see 
the huge difference in depth between the two wells. The inset shows a closer look to the minimum 
located at q = 0, to provide a better idea of what the energy barrier is in this case. 

population, that we have to make it easier for them to cross the barrier. The effect of a 
nonvanishing third moment in the noise term has been previously examined in the case of a 



deterministic chaotic "noise" [17j|. But in this case the weakness of the asymmetry allowed 
to treat it as a perturbation. 

The relaxation time to go is given by a = (L* — 1) , so we have all the ingredients to 
calculate the decay time r, and we can now examine if our system will be resonant to the 
external signal f(t) = A cos (out + <fi). Here is an arbitrary phase and A is the amplitude 
of the signal. In our population dynamics model it is introduced as a periodic modulation 
of the competition — ► + A cos {ojt + 0), and we will assume that the amplitude A of 
this signal is very small so we can consider it as a perturbation. This perturbation implies 
an additional time-dependent term in the Hamiltonian H p = n 3 Acos(um~ 2 t + (p)(p — l)p 2 q 3 , 
in non-dimensional variables. The corresponding action reads 

S p = nA cos(ujn~ 2 t Q ) / cos{uon~ 2 t)(p— l)p 2 q 3 dt— nA sm(un~ 2 t ) / sm(ujn~ 2 t)(p— l)p 2 q 3 dt, 
Jo Jo 

(10) 
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where we have expressed the phase as a function of the signal initial condition = u>n~Ho. 
Now, after assuming a quasistationary variation of the external signal (adiabatic approxi- 
mation), we can reduce the action to 

S p = nAcos(ujn- 2 t) [ (p - l)p 2 q 3 dt, (11) 
Jo 

where we have modified the notation t — > t. To compute action we need to solve 
the classical equations of motion, q = —d p H and p = d q H, in the instanton line Eq. (JZj); 
for the first equation we get 2dt = (^L*q — \/L*q 3 (l + L* — q)^j dq. Plugging this last 
expression together with Eq. (JJJ) (note that now e = and L = L*) into Eq. (|TT|) we obtain, 
for transitions from q = to q = 1 

f 1 y/ L*q 
S p = nA cos(un~ 2 t) / — — = 



l 2(1 + L* - g) 3 / 2 

nAcos(ujrr 2 t) (\ — \/L* ■ arccsc L* + 1 j J = nA\ cos(ujn~ 2 t), (12) 
and for transitions from q = L* to q = 1 

S p = nA cos(^" 2 t) £ 2(1+ ^ g)8/a <*g = 

nA cos(um~ 2 t) ^1 — L* + vT* arctan ^\/Z*J — arccot ^VL*J J = —nA 2 cos(un~ 2 t), (13) 

where the effective amplitudes Ai and A 2 were defined in order to preserve their positiv- 
ity. The ratio between both actions \S P (L* — > 1)/S P (0 — > 1)| ~ 48 shows that the effective 
amplitude of the external signal is very different at the two stable fixed points, being much 
stronger in q = L* than in q = 0. It seems that this mechanism is reminiscent of that of 
noise amplification, that is responsible in turn of the asymmetry of the potential in Fig. (j2J) . 
So we can finally write the time-dependent escape rates 

n(0 -»• 1) = exp[5 + S p (0 -»■ 1)] « r(l + tiAx cos(a;n- 2 t)), (14a) 
r 2 (L* -> 1) = exp[5 + S P (L* -> 1)] « r(l - nA 2 cos^rt)). (14b) 

Now we can treat our problem as a two-state system governed by the master equation 
P ± (t) = -iy T (t)P ± (t) + W / "±(t)P T (t), where P+ (P_) denotes the probability that the system 
occupies the state q = L* (q=0) at time t, and W+- denote the transition probability 
densities, given by 

W- = W(L* -> 1) = — — « r(l + nA 2 cos^n" 2 *)), (15a) 

r(l — Wi4 2 cos(cjn~^t)) 

W + = W(0 -»• 1) = — ; 7 ^ « r(l - nAi cos(cm~ 2 t)), (15b) 



where r = r _1 . From now on, our derivation of the stochastic resonant behaviour will be 
parallel to that of Ref. so we will not do the calculations in detail and we refer the 
interested reader to this reference. We can use the normalization condition P + + P = 1 
together with the transition probability densities Eqs. (|15a|) and (jl5bj) to solve the two-state 
master equation to first order in A. This solution can be used to obtain the system response 
(q(t)\q(to), to) = J xV(x,t\xo,to)dx, where V(x,t\xo,t ) = P + (t)8(x — L*) + P_(t)5(x). In 
the asymptotic limit we find 

i / / m / \ \ L* L* n 3 r r uj / uj \l . . 

hm (q(t )\q{to),t ) = - -(A 1 +A 2 ) , cos —t — arctan ( - - J , 16 

t -*-oo 2 2 V4r 2 n 4 + u 2 \-n A \2m A /\ 

and we can appreciate that the amplitude of the periodic part of the system response under- 
goes a resonance for a finite value of n. To see this more clearly consider the n-dependent 
part of the amplitude of the system response 

n 3 r n 3 (L* — 1) exp(— nR) , , 

A n = , = , (17) 

V4r 2 n 4 + uj 2 ^A(L* - l) 2 exp(-2nR)n 4 + u 2 
where R = Sq/h is a constant independent of n. This function attains its maximum at 
the value of n solving the transcendental equation 4(L* — l) 2 n 4 = exp(2nR)(nR — 3)u 2 . 
In order to get a more graphical idea of the phenomenon we have depicted A n versus n in 
Fig- © for three values of the frequency u = 10~ 3 , 10~ 4 , 10~ 5 . In this figure one can see how 
the resonance becomes stronger and appears for higher values of n as the forcing frequency 
decreases. 

Let us now point out three final and very important remarks. The first is that the 
resonance is present only if the external signal is introduced modulating K4, K 3 , or K 2 . 
A simple analysis reveals that modulating K~i, K 2 , K% shifts the cubic dependence in n in 
A n to ra°,n 1 ,n 2 respectively, so the possibility of resonance is lost (at least, if we do not 
introduce n-dependence in the external signal) for K\. Also, this is not the unique example 
of stochastic resonance that might appear in a chemical system: different parameter values 
and reaction sets will provide new examples of this phenomenon. The assumptions made 
in this work were introduced to facilitate the analytical assessment of the problem. Finally, 
the improvement in the resonance observed in Fig. (JHJ) for decreasing frequencies is not 
unbounded; if we assume too large values of n then we could not consider the action (JTTJ) as 
a perturbation, and so the linear response theory would no longer be valid. 

In summary, we have shown that it is possible to find stochastic resonance in reaction 
kinetics purely due to the presence of intrinsic noise generated by the discrete character of the 
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FIG. 3: Nonmonotonic behaviour of A n as a function of n for three different values of the forcing 
frequency: u = 10 -3 (solid line), ui = 10~ 4 (dashed line), and u> = 10~ 5 (dotted line). We can see 
how the resonance appears for higher values of n and becomes stronger as the forcing frequency 
decreases. 

reactants. We have presented our model as a population dynamics problem, a type of system 
that is usually subject to external periodic forces, as seasonal variation. A possible resonant 
coupling between phenotype selection in a biological species and periodic environmental 
evolution has been suggested recently Q|, so it would be very interesting to know if not only 
the phenotype but also the population itself can undergo some sort of stochastic resonance, 
and in particular the one reported here. This could have a serious impact on the possibility 
of extinction of a population, since, as we have seen, internal fluctuations can drive a system 
from a state with a large population to an empty state. 

This work has been partially supported by the Ministerio de Education y Ciencia (Spain) 
through Projects No. EX2005-0976 and FIS2005-01729. 
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